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Abstract In calculations of the inspiral of binary black holes an intermediate approximation is 
needed that can bridge the post-Newtonian methods of the early inspiral and the numerical relativity 
computations of the final plunge. We describe here the periodic standing wave approximation; A 
numerical solution is found to the problem of a periodic rotating binary with helically symmetric 
standing wave fields, and from this solution an approximation is extracted for the physically relevant 
problem of inspiral with outgoing waves. The approximation underlying this approach has been 
recently confirmed with innovative numerical methods applied to nonlinear model problems. 



The computational study of the inspiral of binary black 
holes is important for the understanding of gravitational 
wave signals, and is of inherent interest as a question in 
general relativity that can be answered only with com- 
putation. It has therefore become the focus of super- 
computer codes that evolve Einstein's field equations for- 
ward in time from initial conditions chosen to represent 
a starting configuration of the inspiralling objects. The 
evolution codes, however, typically become unstable on a 
timescale (set by the size of the hole) short compared to 
a full orbit. Reliable calculations of the final plunge are 
now feasible, the merger and ringdown of the final black 
hole fate of the system are handled well with perturba- 
tion theory 1], and the early inspiral is well approximated 
with post-Newtonian computations What cannot be 
handled well is the intermediate phase of the inspiral, 
that late epoch during which nonlinear effects are too 
strong for the post-Newtonian approximation, and too 
many orbits remain for numerical relativity to be stable. 

It has long been recognized that the basis of an ap- 
proximation scheme should be the slow rate of inspiral, 
the small ratio of the orbital time to the radiation damp- 
ing timej^, 0|- Through an adiabatic treatment of the 
slow inspiral, such an approximation could give answers 
about the radiation and rate of inspiral in the interme- 
diate epoch. And when the rate of inspiral became too 
rapid, the intermediate approximation could hand the 
problem off to numerical evolution codes to do the fi- 
nal orbit and plunge, and could supply the ideal initial 
data to those evolution codes. The need and the concept 
for an intermediate approximation have been clear, but 
such an approximation has not been easy to implement. 
Along with several colleagues 0,0. 013 we have based an 
approximation of slow inspiral on a numerical computa- 
tion of no inspiral. That is, we seek a numerical solution 



of Einstein's equations for binary objects that are in cir- 
cular periodic motion, and whose "helically symmetric" 
fields rotate rigidly with the source objects. 

The universality of gravitation suggests that the un- 
changing motion of such a system is not compatible with 
outgoing radiation, and this intuitive suggestion is con- 
firmed by the mathematics of the theory. We, there- 
fore, seek a helically symmetric solution for the sources 
coupled to standing waves. In a linear theory standing 
waves would be a superposition of half-ingoing and half- 
outgoing solutions. In linear theory, one could (though 
without motivation) solve the standing wave problem. 
From the fact that solution is half the superposition of 
the ingoing and outgoing solutions, and from the rela- 
tionship of the ingoing and outgoing solutions, one could 
extract the outgoing solution. The crux of our periodic 
standing wave method is that even for highly nonlinear 
binary inspiral fields there is an "effective linearity." The 
standing wave solution, to good accuracy, is half the sum 
of the outgoing plus ingoing solutions despite the non- 
linearities. In general relativity, therefore, we should be 
able to solve the standing wave numerical problem and 
extract an approximation to the outgoing solution. 

It is important to understand why effective linearity 
can be correct for inspiral. In the strong-field regions 
very close to the sources, the solution is very insensitive 
to the distant radiative boundary conditions (ingoing, 
outgoing, standing wave). In this near-source region a 
superposition of half the ingoing and half the outgoing 
solution gives a good approximation solution, because 
it amounts to averaging two samples of the same thing. 
In the wave zone where the outgoing and the ingoing 
solutions are very different, the fields are weak enough 
that nonlinear effects are negligible, and once again we 
can superpose. The separation of the strong-field region 
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from the boundary-influenced region should be a clear 
separation unless the sources are rotating very close to 
c, in which case the wave zone will start just outside the 
sources. It is, however, not expected that ultrarelativistic 
source motion can occur during the slow inspiral epoch 
of motion. 

We have recently been able to confirm effective linear- 
ity. (Technical details are given in This confirma- 
tion has been achieved with a model problem, since the 
validity of effective linearity can only be carried out in 
a model problem. In general relativity, there will be no 
"true outgoing" solution available for confirmation until 
numerical evolution codes are fully developed. In addi- 
tion, the numerical features of the helically symmetric 
standing wave calculation pose new challenges very dif- 
ferent from those of evolution codes, and are best resolved 
in the simplest context possible. 

Our model problem is a nonlinear scalar field coupled 
to point-like sources in Minkowski space, and satisfying 
the field equation 

^'•a /sg"^ + AF = V^* - 4^?* + F = Source . (1) 

Here the source is taken to be two points of unit scalar 
charge in orbit around each other at angular frequency Q, 
and at radius a. The velocity parameter for the system 
/3 = afl /c is taken to be of order unity, representing the 
strong-field tight binary for which post-Newtonian ap- 
proximations are inadequate. Our methods, however, re- 
quire that P not be too close to unity. More explicitly, our 
method works best if the radiation is quadrupole domi- 
nated, and our calculations focus on values of (3 from 0.3 
to 0.5. We expect that the approximation of slow inspi- 
ral will break down before this assumption breaks down, 
and numerical relativity evolutions will take over the job 
of tracking the last part of an orbit and the subsequent 
plunge and merger. 

The term F contains the nonlinearity in our model 
theory, and we have found the following form, with pa- 
rameters A and ^Pq, to be very useful: 
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A crucial feature of F is that like the nonlinearities of 
general relativity, it is very large near the sources, and 
becomes negligible far from the sources. The A multiplier 
allows us to vary the strength of the nonlinear term, and 
the ^'o parameter allows us to vary the profile of the 
nonlinearity in the strong field region. 

Our problem is defined by Eqs. |^ and (0), and the 
source motion at angular frequency f2 in the equatorial 
plane. As described in spherical coordinates, helical sym- 
metry can be imposed on the solution r, 9, </>) by re- 
stricting to solutions of the form ^'(r, 9, ip), where ip is the 
comoving azimuthal coordinate (p—VLt. By restricting the 



solution in this way, we have eliminated the possibility 
of "evolution." For such helically symmetric solutions a 
change in time by At is the same as a change the az- 
imuthal angle A(/) = — f2At. With this suppression of 
evolution, we have eliminated the sorts of instabilities 
that develop in evolution codes. But we have introduced 
new difficulties. 

These new difficulties can be seen immediately in the 
form of the helically restricted nonlinear scalar equation 
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= Source -F(*) = ct(^'). (3) 



The principal part of this quasilinear equation is 
"mixed," elliptic inside a cylinder at rsin6' = c/fi, and 
hyperbolic outside that cylinder. The problem is to be 
solved with radiative conditions (ingoing, outgoing, or 
standing wave as described below) on a spherical surface 
at large distances. Well posed problems in physics typi- 
cally supply cauchy data on open surfaces to hyperbolic 
equations, and Dirichlet or Neumann data on closed sur- 
faces to elliptic equations. Our model leads to a bound- 
ary value problem with "radiation" conditions on a closed 
surface surrounding a mixed problem. Though unusual, 
our problem is intuitively well posed, and passes a com- 
putational test: we have found no fundamental difficulty 
in solving models of this type numerically. Furthermore, 
a careful analysis^^ of a closely related problem proves 
that solutions exist and are stable. 

"Standing wave" solutions - half ingoing and half out- 
going -are at the heart of our method, but there is not 
an obvious definition of standing wave solutions in a non- 
linear theory. Our procedure is to find the outgoing 
and ingoing Green functions for Eq. In princi- 
ple, we can then iterate to find a solution of Eq. Q. The 
iteration 
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if it converges, gives V&out, our nonlinear outgoing solu- 
tion (and similarly for 4"^), while the convergent result 
of 
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is what we mean by our nonlinear standing wave solution, 
"^std- The standing wave solution '^std is fundamentally 
different from {'^out + ^in)/2, but if effective linearity is 
correct, the two are very nearly equal. (Note: In practice, 
for strong nonlinearities, the direct iteration described 
above must be replaced by Newton-Raphson iteration.) 

We have previously!^ solved the model problem of 
Eq. lO with more-or-less straightforward finite differ- 
encing and direct matrix inversion. (The mixed nature 
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of the partial differential equations prevents the use of 
such efficient techniques as explicit relaxaton.) This ap- 
proach was successful (iterations converged) for models 
with a limited range of source velocities and nonlinear- 
ities. More recently we have developed an innovative 
numerical method that gives remarkably good results, 
with very little computational cost, and that might be 
useful in problems other than ours. Our new method is 
based on three elements. First, we use "adapted coordi- 
nates," comoving coordinates that conform to the geom- 
etry of the problem. Near the source points our coordi- 
nate surfaces approach those of source-centered spherical 
harmonics; far from the sources the coordinates become 
spherical polar coordinates centered on the midpoint of 
the orbit. Our specific choice of adapted coordinates is 
two-centered bipolar coordinates Xi^^^ (pictured in the 
equatorial plane of the orbit in Fig.^, which asymptoti- 
cally approach spherical coordinates r, 6, </). As discussed 
in this choice is not the most computationally ef- 
ficient possibility, but it has the advantage of relative 
simplicity. 




mathematics evaluated on the discrete numerical grid. 
We have found that this approach does not work at all 
well for our purposes. Because the monopole is so much 
greater than the quadrupole in the wave zone, project- 
ing out the quadrupole component with the continuum 
multipole gives a large error. We have therefore used the 
multipoles that seem ideally suited to decomposition on 
the angular grid with ne x grid vertices. We view the 
values of our solution as vectors ^'(6^, 4>j) in a space of 
dimension uq x At large x the angular Laplacian 

^ _2 d_ f . _d_\ 1 

- sin e ae y'"" ae)^ sm^ e a$2 

can be implemented as a finite difference operator, on 
the HQ X n$ space, that is self-adjoint with respect to 
the finite difference equivalent of integration over solid 
angle0. The eigenvectors of this self-adjoint operator 
are approximately Yim{&i,^j), but the eigenvectors are 
exactly orthogonal so that the projection of the radiative 
multipoles is not affected by the large monopole. 
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FIG. 1: TCBC coordinates in the equatorial plane. 

The second element of our method is to expand in mul- 
tipoles (spherical harmonics) of the adapted coordinates. 
If the separation between the sources is much larger than 
the size of the source, then the only feature of the source 
that is important is the "local" monopole (monopole in 
source-centered coordinates). In the wave zone, on the 
other hand, only the quadrupole of the field matters for 
the radiation (unless the source motion is highly rela- 
tivistic). This suggests that a multipole expansion in the 
adapted coordinates need retain only the monopole and 
quadrupole. This turns out to be true for source speeds 
< 0.3 c. For somewhat larger speeds, good accuracy re- 
quires that the hexadecapole be kept in addition. This 
severe filtering of the multipoles reduces the computa- 
tional burden of solving the problem, but more impor- 
tant, it eliminates the numerical noise at short angular 
scales that we found for straightforward finite difference 
computations. 

An "eigenspectral" treatment of multipoles is the third 
element of our method that is innovative. A straightfor- 
ward approach to dealing with multipoles would be to 
use y?m(Oi, $j) on an angular grid of the O and <& coor- 
dinates, that is, to use the multipoles of the continuum 



FIG. 2: A comparison of exact (series) linear outgoing so- 
lutions with eigenspectral solutions. The radiation field (^P 
minus its monopole) is plotted against the comoving coordi- 
nate Z, the distance from the center along a line through the 
source points. Eigenvectors were found on a 40 x 80 grid for 
one quadrant, and 16001 radial grid points were used, with an 
inner boundary condition at 0.02 a, and a Sommerfeld outer 
condition at 80 a. 

Questions about the validity of our eigenspectral 
method apply just as much to the linear version as to 
the nonlinear version of our model problem, and for the 
linear problem exact (infinite series) solutions exist for 
comparison. Figure |21 shows such comparisons for mod- 
els with source velocity aVl = 0.4c and aVl = 0.5c and 
shows multipole filtering that allows either the monopole 
plus quadrupole {£max — 3), or with the hexadecapole 
included {imax = 5). As the figure indicates, and as 
should be expected, the number of multipoles that need 
be included to achieve a given accuracy increases with 
increasing source speed. The minimal (max — 3 mul- 
tipole set may be adequate for approximate results at 
aH. — 0.4c and gives excellent results for aH. — 0.3c (not 
shown). There is, of course, no need to limit the num- 
ber of multipoles retained to such a small set, but as the 
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number of multipoles approaches the maximum number 
of eigenvectors that can be found for the angular grid, 
the smoothing effect of multipole filtering is lost and nu- 
merical noise becomes significant. A large number of 
multipoles would be necessary only for source speeds af2 
close to c. As already mentioned, it is unlikely that such 
motion will satisfy the slow inspiral condition for which 
our approximation is designed. 

Other confirmations of the validity of the method, es- 
pecially for nonlinear models, have been carried out, and 
are reported in Ref. ■ We focus here on the most im- 
portant question that can be answered with these models 
and numerical methods: Does effective linearity work? 
Can we extract a good approximation to the outgoing 
nonlinear problem from the sort of standing wave compu- 
tation we will be limited to when dealing with Einstein's 
theory? Figure gives strong evidence that we can. For 
the chosen parameters A = —15 and \E'o — 0.15, non- 
linearities are significant, strong enough to reduce field 
strength by around two-thirds. The outgoing and stand- 
ing wave solutions were each computed by the Newton- 
Raphson version of the iteration in Eqs. Q,©. 

An outgoing approximation is extracted from the 
standing wave solution by the following steps: (i) In the 
outer region the solution is matched to a general linear 
standing wave solution of half-ingoing and half-outgoing 
waves; an outgoing wave is extracted as if the problem 
were linear, (ii) In the strong-field region very close to 
the source, the standing wave solution itself is taken to 
be the outgoing approximation, (iii) The two solutions 
are blended over a narrow intermediate range of radii. 

In Fig. 13 the computed outgoing nonlinear solution is 
shown as a solid curve. The data-type points representing 
the outgoing wave extracted from the standing wave so- 
lution show how good the approximation is. We have run 
models with much stronger nonlinearity and have found 
equally good, or better, agreement of the true outgoing 
solution and the extracted approximation. The validity 
of effective linearity should, in fact, become questionable 
not for stronger nonlinearity, but only for physically im- 
plausible high source velocity. 

In addition to confirming effective linearity, computa- 
tion with the model has also allowed some early insights 
about sensitivity to source details. By varying the multi- 
pole content of the inner boundary data we explored the 
impact of source structure on the radiation field. The 
result (detailed in Ref. ^|) is in perfect accord with 
physical intuition; the radiation is insensitive to source 
structure unless the source size becomes comparable to 
the separation of the sources (i.e. , unless the moments 
ascribable to the structure of the individual sources are 
comparable to the quadrupole moment due to the sepa- 
ration of the mass points). The equivalent question for 
Einstein's theory is more difficult, but we should be able 
to give clear quantitative answers. 

The next steps in our program start with linearized 



gravity in the harmonic gauge. We have already done this 
with a finite difference code; work on applying the eigen- 
spectral method is underway, and no significant prob- 
lems are anticipated. The infrastructure of the linearized 
gravity will provide much of what is needed for post- 
Minkowskian computations, since the structure of the lin- 
ear operators (the analogs of the £s in Eqs. lEl)-©) will 
be the same as for linearized gravity. The final step to the 
full Einstein theory will, in the same way, be based on the 
numerical infrastructure developed for post-Minkowsian 
models. 
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FIG. 3: The computed nonlinear outgoing solution compared 
with an approximate outgoing solution extracted from the 
computed nonlinear standing wave solutions. Grid parame- 
ters are the same as for Fig. |21 and ^max = 5 was used. 
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